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The expression "spin glasses" was originally coined to describe metallic alloys 
of a non-magnetic metal with few, randomly substituted, magnetic impurities. 
Experimental evidences where obtained for a low temperature "spin-glass" 
phase characterized by a non-periodic freezing of the magnetic moments (the 
spins) with a very slow, and strongly history dependent, response to exter- 
nal perturbations (this later aspect lead more recently to many fascinating 
developments). The theoretical analysis of this phenomenon lead to the cel- 
ebrated Edwards-Anderson model [1] of spin-glasses: classical spins on the 
sites of a regular lattice with random interactions between nearest neighbor 
spins. However, after more than thirty years of intense studies, the very na- 
ture of the low temperature phase of the Edwards-Anderson model in three 
dimensions is still debated, even in the simple case of Ising spins. Two main 
competing theories exist: the mean field approach originating from the work 
of Sherrington and Kirkpatrick [2] , and the so-called "droplet" or scaling 
theory of spin glasses. 

The mean field approach is the application to this problem of the conven- 
tional approach to phase transitions in statistical physics: one first builds a 
mean field theory after identifying the proper order parameter, solve it (usu- 
ally a straightforward task) and then study the fluctuations around the mean 
field solution. Usually, fluctuations turn out to have mild effects for space 
dimensions above the so-called upper critical dimension (up to infinite space 
dimension, where mean field is exact). Below the upper critical dimension 
fluctuations have major effects and non-perturbative techniques are needed 
to handle them. The second item of this agenda (solving the mean field equa- 
tions) led, with spin glasses, to severe unexpected difficulties, and revealed a 
variety of new fascinating phenomena. The last step is the subject of the so 
called "replica field theory" , which is still facing formidable difficulties. 

These notes arc an introduction to the physics of the infinite range ver- 
sion of the Edwards-Anderson model, the so-called Sherrington-Kirkpatrick 
model, namely a model of classical spins that are not embedded in Euclidean 
space, with all pairs of spins interacting with a random interaction. If there is 
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no more debate whether Parisi famous solution of the Sherrington-Kirkpatrick 
model in the infinite volume limit is correct, much less is known, as mentioned 
before, about the Edwards- Anderson model in three dimensions, with numer- 
ical simulations as one of our main sources of knowledge. It is accordingly 
important to test the various methods of analysis proposed for the Edwards- 
Anderson model, in the Sherrington-Kirkpatrick model case first. 

In a first part, I motivate and introduce the Edwards- Anderson and 
Sherrington-Kirkpatrick models. In the second part, I sketch the analytical 
solution of the Sherrington-Kirkpatrick model, following Parisi [H O [SI El- 
I next give the physical interpretation of this solution. This is a vast subject, 
and I concentrate on the major points and give references for more develop- 
ments. The third part presents the numerical simulation approach and com- 
pares some numerical results to theoretical expectations. The last part, more 
detailed, is about the specific problem of finite size effects for the free energy, 
which is interesting for both theoretical and practical point of views. I have 
left aside several very interesting aspects, like the problem of chaos [H [9] the 
TAP approach [10] (see [11] for numerical results) and the computation of the 
complexity [12] . 

There are many books and review articles about spin glasses and related 
phenomena: One may start with the text book by Fischer and Hertz |13j . 
the review by Binder and Young [14| . which gives a very complete account 
of the situation in 1986 both experimental, theoretical and numerical with 
many detailed analytical computations, and (at a higher level) the book by 
Mezard, Parisi, and Virasoro [15]. More recent references include [16] and [17] . 
The recent book by dc Dominicis and Giardina [TH] gives a very compete 
exposition of the replica field theory. Reviews on various aspects of the physics 
of spin glasses can be found in [lH [20l [21]. The non equilibrium behavior 
of spin glasses have been the subject of intense work during recent years, 
see [221 [231 [24] for reviews. 

1 Introduction 

The Edwards-Anderson Ising (EAI) model is the classical Ising model with 
quenched random interactions. The Hamiltonian is 

= - ^ Jtj(Ti<7j - H^ai , (1) 

<ij> i 

where the variables ct^ = ±1 are Ising spins living on the sites of a regular 
square lattice in d dimensions. H is the magnetic field. The interaction involves 
all nearest neighbor pairs of spins, with a strength Jj^^ that depends on the 
particular link < i,j >. The Jij^s are quenched variables, namely they do not 
fluctuate. They have been drawn (once for ever) independently from a unique 
probability distribution P{J) with mean Jq and square deviation J^. In what 
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follows, I will consider the case Jq — and choose = 1. The later choice is 
just the choice of the unit of temperature. 

For a fixed set of Jij 's, denoted by J^, one can compute thermodynamic 
average values, for example the average internal energy Ej =< E >, where 
the symbol < • • • > denotes the thermodynamics average (with given J') . In 
general the results depend on J', and should be averaged over the J' distri- 
bution. I note by this disorder average. 

This model was proposed [1] as a simple model that captures the essen- 
tial of the physics of spin glasses, and in particular of magnetic spin glasses. 
Those are alloys with magnetic impurities randomly distributed inside a non 
magnetic matrix. The interaction between the impurities is the RKKY inter- 
action that oscillates with the distance. Since the positions of the impurities 
are random, their interactions are random too and one is led to the Edwards- 
Anderson model (see e.g. (TJ [T31 E] for more than this abrupt summary). 
If this very simple model is to explain convincingly the behavior of real spin 
glasses, its physics should obviously not depend (too much) on the specific dis- 
order distribution used, as soon as it has zero mean and unit square deviation. 
This disorder distribution is usually Gaussian (this leads to simpler analyt- 
ical computations) or binary with values ±1 (this leads to faster computer 
programs, using the multi-spin coding technique). 

The same model can be generalized to vector spins with two or three com- 
ponents. Those are called the XY and Heisenberg Edwards Anderson mod- 
els respectively. Most real spin glasses are indeed Heisenberg spin glass with 
anisotropic interactions (due to the lattice structure, the interaction is not the 
rotational invariant cTiCTj). 

In the zero magnetic field case {H = 0), there is now an agreement (based 
mostly on numerical simulations) between spin glass physicists that the EAI 
model has no transition at finite T in two dimensions, and a transition for three 
dimensions and above 0. In the low temperature phase the spins are frozen, 
with < Gi >^ 0, but with a random pattern, and < ctj > is of order 
l/y/N , where N is the number of spins. There is accordingly no spontaneous 
magnetization. There is no hidden magnetization either: the spins follow no 
periodic pattern (like e.g. in an anti-ferromagnet), and it is furthermore very 
likely that the spin orientations arc completely reshuffled as soon as one varies 
the temperature H. This is summarized by the statement: The EAI model is 
not a disguised ferromagnet. 

The statement < Ui >^ needs to be made more precise. For a finite 
system, it means that < Ui > is nonzero when observed over a time scale 
t << te.rg{N), where terg{N) is called the ergodic time. At some later time 
t, the system will eventually tunnel from the current equilibrium state, to a 
state where < ct; > is reversed. However the dynamics of the EAI model is 
extremely slow and terg{N) is enormous in the low T phase, as soon as one 



^ Tlie value of the lower critical dimension, between 2 and 3, is not yet settled. 
^ This is the so called temperature chaos effect. 
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has few hundred spins, to the point that with any real spin glass, this time is 
much larger than the duration of any experiment. 

In absence of magnetization, we have to find another order parameter for 
the glass order, one uses the so called Edwards- Anderson parameter, namely 

= hm lim -^-S^ < ai{tn)ai{t + t^) > , (2) 

i 

that is nonzero below the spin glass transition and zero above. On a finite 
system, C{t) = l/N"^- < ai{to)ai{t + to) >, will rapidly decrease from the 
starting value C(l) = 1, have a plateau of height qea for a long time (of order 
terg) bcforc dropping. As A'' grows, the plateau becomes longer and longer. 
The so defined qea turns out to depend weakly on J for large systems. This 
is a so called self averaging quantity (see Subsec [2?3)) . 

One may also consider two independent copies of the system, {cTj-^''} and 

(2) 

{(tI '}, with the same disorder instance J' (such copies are called real replica, 
or sometimes clones, in order to distinguish them from the replica of the 
replica method, see Subsec. I^T^ . and consider the probability distribution 

P^(g)=<%-^E'^'"^f')>' (3) 

i 

namely the probability distribution of the overlap q = 1/-/V crp'o-p' be- 
tween the two systems. A nonzero qea corresponds to two peaks in Pj{q) 
centered at q = ^Qea- We will see later that, at least for the SK model, 
Pj{q) has more structure than this double peak. 

Clearly, qea =< M >^ for a ferromagnet in d dimensions (still at zero 
magnetic field), where < A/ > is the spontaneous magnetization. On a finite 
system, P{q) is made of two Gaussian centered around ± < M >^. Between 
the two peaks, P{q) is exponentially small [25], of order exp(— 2y^L'^~^), where 
A is the interface tension, and L the linear dimension of the system. 

In ferromagnets, the transition happens when the magnetization M ac- 
quires a nonzero expectation value, and accordingly the susceptibility x = 
A'^ < > becomes infinite as T decreases towards Tg. In spin glasses the 
spin susceptibility stays finite as T — > Tc, but the spin-glass susceptibility 
XSG = N < q^ > diverges. There is however a difference: in Ising ferromag- 
nets, one can define a connected susceptibility A^(< Af ^ > — < M >^) 
that is finite in the low T phase. For Ising spin glasses the analog connected 
susceptibility diverges, as A^ ^ oo in the whole low T phase 



In this formula < Af > is to be interpreted as an average restricted over a time 
interval of length << terg{N), or simply as < \M\ >. 
* It is usually defined as "^^i j (< '^i'^j > — < o^i >< aj >y, sometimes simply as 
N{< g2 > _^g^2-j ^ 
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2 The Sherrington— Kirkpatrick model 

This model [2] is a simplified version of the Edwards-Anderson model, where 
all spins are directly coupled. The Hamiltonian is 

^ = - E 7=^^'^.-^E'^- (4) 

l<i<j<N ^ i 

where N is the number of spins. The Jij 's are again drawn from a unique 
probability distribution P(J') with zero mean and unit square deviation. The 
factor will ensure a finite limit for the internal energy per spin eN{T) 

as iV ^ oo. It is in some sense a model in infinite dimension since a given spin 
is coupled to N spins, with A'^ — > oo in the thermodynamic limit. The model 
has "infinite connectivity" . The infinite connectivity will ensure that the mean 
field method gives the exact result. As we will show, the SK model can be 
solved exactly in the thermodynamic limit, and the result is independent of 
the choice made for the disorder distribution P{J)- 

The standard method to solve this model is the famous replica trick 0. 
One is interested in computing the average free energy —(iFj = InZj-, where 
(3 = 1/T is the inverse temperature, in units such that fcs = 1. This is 
done by considering the partition function of n identical un-coupled copies of 
the system, with the same instance of the disorder J . Such copies are called 
replica. The partition function is simply the nth power ol Zj, Zj-. Continuing 
this function defined for integer n > 1 to a function of the real variable n, one 
has, at least for finite N , 



_ 2 

pF = \-siZj= lim . (5) 

n— »o n 



Z^j can be written as 



Z-j^Tr^-^eM^m^j^)) (6) 
TrN exp(-/3(7^j(a(i)) + • • • , 

where a^^^ represents the A^ spins of the first replica (namely {ct'^"*}), ct'-^-' 

[nl 

the A^ spins of the second replica, . . . , and Tr\j is the trace over the nN spin 
variables. The average over the 71(71-- l)/2 Jij variables are independent, and 
the disorder average factorizes as a product of terms of the form 



* There is an alternative method called the cavity method, see [151 126] . 
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where 

n 

Y _ ^(1)^(1) I ^(2)^(2) I ,Jn) {n)_^ {a) {a) 

-^i,3-'^i ^] ^ ~ Z^'^i '^j ' l^) 

a=l 

and the [J]p's are the successive cumulants of the disorder distribution P{J), 
[J], = J (9) 



[jL = {j^jy 



[j]^ = {j^jY-i{j-jY 

Since we have J = and J"^ = 1, one is led to the result 



The neglected terms in the exponent are of order 1 /N'^ and will not con- 
tribute to the thermodynamics limit and accordingly the physics is inde- 
pendent of the disorder distribution. One obtains for the disorder averaged 
partition function: 



i<j fc=i i b=i 

= ^^Klv + /5^E E ^ 



/?2 

The average over the disorder has been performed analytically, but now 
the n replicas are coupled. In order to proceed further, one uses the formula 



j dq exp( ^) exp(g/32x) = cxp(^) , (12) 

introducing n(n — l)/2 auxiliary real variables qa.t {o, < b), with the result 



They are altogether absent if the distribution is Gaussian. 



Infinite range models of spin glasses 



7 



a<b a<ib 

+ E E -.^"^-f ^ +^hy^T. + ("^ - -')t) ■ 

The variables ga.h have been defined for a < 5. In the foUowing. it wiU be 
sometimes convenient to define qa,b for a> b also, as qa,b = qt.a and ga,a = 0. 
For a given replica index a, the trace over the spins factorizes as 



TV 



/3'E'^-E-l"^-f^+'3^EE-. 



(b) 



J|Tr^(a) e °<f' * * ''=1 (14) 

n 

where 5^"' is any of the a^'^'^^s, or alternatively the spin of a single spin system, 
which is replicated n times. Thus0 



a<b \ a<b 

n 

N log (Tr [exp(/32 ^ g.^.^t-) ^t^) + ^ S^"' )]^ 



a<fc b=l / 

The formula has a particularly simple form, with all dependence in N explicit 



^J=[I[\I^J '^'^-^l cM-Nf3A{{qa,b})) . (16) 

a<b 

Up to now our derivation is exact (for a Gaussian disorder) in the n ^ 
limit. We now make the saddle point (or steepest descent) approximation, 
which gives the correct N oo behavior. Assuming the existence of a unique 
absolute maximum of the integrand at location {^f f"}; '^'^'^ ^^^^ 

Z^^cM-NPA{{q!^})). (17) 



From now on, we omit the terms in the exponent, since they do not contribute 
in the n — > limit. 
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By assumption, all partial derivatives of yl, dA/dqa^, are zero at the 
saddle point, and the matrix of the second derivatives, the Hessian, has only 
nonnegative eigenvalues. The free energy of the original SK model is thus 
simply related to the saddle point of the qa^b integral representation of the 
partition function of the n times replicated model (with a reckless interchange 
of limits), 

/= lim 1^ = lim i^aqf ri) . (18) 

Note that the saddle point equations can be written as self consistent 
equations, involving a single spin system 



Tr5[5(-)5('') exp(/?2 + (3HY,S^'^)] 
.ir = ^ ■ (19) 

a<b b=l 

The replica method is not limited to the evaluation of the free energy. It 
can be used to compute the average (thermodynamic average and disorder 
average) of any function of the spins. Let 0{a) be a function of the spins 
{ci}, we have by definition, for any disorder sample: 

<0(,)>.^'-""'°"''"-'-^'""». (20) 

Multiplying both numerator and denominator by and letting[f| n — > 

0, we have 



< 



0{a) lim Tri"l0(a(i))exp(-/3(H'"1)) . (21) 



The right hand side is of a form whose disorder average is readily computed 
using (fTO| . The method is extended readily to the disorder average of products 
like < 0{(j) >< V{<j) >. Introducing two real replica, we have indeed 



Tri'l0(a(i))7'(a(2)^) exp(-/3(H[^l)) 



<0{cj)><V{a)>^ ^ ' • ^ J" , (22) 

'J 

Multiplying both numerator and denominator by and letting n — > 0, 
we have 



< 0{a) >< V{<j) lim rr[,"lC'(cr(i')7'(cr(2))cxp(-/3(7^["l)) . (23) 
which is again of a form whose disorder average is readily computed using ()10p . 



* We did add n ~ 1 replica to the one real replica and then let n ^ 0. 
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2.1 Interpretation of the qa,,b variables 

The Qa^b variables have been introduced formally. Their value at the saddle 
point has a simple interpretation [57] in terms of the overlaps between real 
replicas as defined in ([3]). Consider two clones cr^^^ and ct'^^ and compute the 
generating function G(j/j2|. 

G(y) = <exp(^^.«af^)>. (24) 

i 

We have, 




Multiplying the numerator and denominator by ^, and letting n 0, 
one obtains after averaging over the disorder 



G(y) = [n / dqa,t] exp - ^ E li>^ (26) 

a<b \ a<b 

+ iVlog(rr|?l[cxp(/32^g,,b5(")5W+/3i/E^*''')]) 



a<h 6=1 



4 



with = 1 since n — > 0. 

At leading order in 1/N, the saddle point is y independent, and accordingly 



< exp(^E-l'^-r') > - -MvP^fS) ■ (27) 

i 

If the saddle point is not symmetric, namely if the q^^'s are not all equal, 
the solution of the saddle point equations is degenerate, and one must aver- 
agf°l over all degenerate solutions, or alternatively all permutations of a and 
b. The correct formula is then 



< -P(^ E -r^-f ^) > = E exp(y/3^'zfl ) , (28) 

i ^ ' a<b 

where n{n — l)/2 is the number of QaM variables. This implies the following 
direct relation between the disorder averaged P{q) and the value of qa,b at the 
saddle point 



^ This is a convenient way to evaluate at once the expression 
< (1/^ Ei erf 'c^f V > for all values of k. 

Since the saddle point solution is degenerate, one must sum over the solutions in 
both numerators and denominators of (I26II. 
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P{q) = Picl)J^<Si<l-T7T.''^ 




(29) 



2 



E'5('?-<r). 



n(n — 1) 



a<b 



The distribution Pj{q), namely before disorder average, has an interpreta- 
tion in terms of pure states. In the rigorous formulation of statistical physics 
(without disorder), pure states give a precise definition of thermodynamic 
phases directly in the infinite volume limit (see [H] for a discussion from a 
physicist point of view). Pure states have the following properties: 

• The clustering property: Inside a pure state, spin correlation functions 
factorize when the distance goes to infinity. For example if a is a pure 
state, one has < a^Cy >a^< cr^ >«< (^y >q, as \x — y\ ^ co. 

• Every translationally invariant state is a convex linear combination of pure 
states. This means that for any A, one has < A Wa < A >q., where 
the weights Wa > 0, with Wa ~ 1, arc independent of A. 

In a ferromagnet at temperature T < Tc, there are two Tdependent pure 
states corresponding to states with positive and negative magnetizations re- 
spectively, that we can call the and the "-" states. A general transla- 
tionally invariant state is a linear combination of these two pure states. For 
example, doing the infinite volume limit (at zero magnetic field) by considering 
a set of finite systems of increasing sizes with periodic boundary conditions, 
one has uj+ = cu- = 1/2. 

Let us assume (departing boldly from rigor) that the same decomposition 
holds for finite systems in the SK model for every disorder sample, namely 
that < A >= Wa < A >a, for any A, where the weights and the states 
are (disorder) sample dependent. Wc introduce Q, the overlap between two 
pure states. 



In a ferromagnet Q+,+ = Q-- =< M >^, and Q+,_ = Q-,+ = — < 
M We now proceed to show that the overlap between pure states is related 
to the clone overlap introduced in Wc do this by considering successive 
moments of q. For the first moment, we have 




(30) 




(31) 
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since cr^^-' and cr'^^) are independent systems. Introducing the pure states, this 
gives 

Q,/3 

For the second moment we have 

77^ EE^"^/5 < 'rf' >„< df^ >„< >p< af^ >p , 

i,j a.,l3 

where we assumed that, since the states a and (3 are pure states, one has 
< aiffj >Q«< ai >a< cFj >a (We pretend that all points are far apart. This 
is clearly not correct for i = j but the error is negligible for large ) . Finally 

<q^ > ^Y.^'^'^Ap ■ (34) 
In general one has, for any integer r, 

< g'' > » E^"^/32a,0 , (35) 

a,l3 

namely (still with disorder dependent Wq's and Qa,i3^s) 

< '^('?- ^E'^^^'^'^^^'^) > = E^"^/"^('?- • (36) 

i a, 13 

This is a remarkable relation between a quantity relative to pure states 
and a quantity that is directly accessible, e.g. with Monte Carlo simulation. 
This relation is quite useful since there is no simple way to characterize pure 
states for spin glasses, in contrast to Ising ferromagnets where a pure state 
can be selected by simply applying an infinitesimal constant magnetic field of 
suitable sign. 
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2.2 Solution of the model 

We look for the absolute minimum of A{{qa,b}), namely the lowest stable 
solution of dA/dqa,b = for all Qa.b with generic n. The solution is then to be 
analytically continued to the limit n 0. One proceeds heuristically by first 
making an ansatz for the matrix qa,6 and generic n, involving a few parameters 
Xi and Qi. Very surprising at first sight, the correct saddle point solution is a 
maximum of A with respect to the Xi's and g^'s (and not a minimum), and 
if several maxima are found, one should usually take the largest . Once a 
candidate solution is found, one should check that this solution is stable (that 
the Hessian has only nonnegative eigenvalues). If no satisfactory solution is 
found, one try a more general ansatz (This is indeed a heuristic procedure). 

The simplest ansatz is the replica symmetric (RS) ansatz [5], namely all 
9a, 6 = 90j and accordingly P{q) = S{q— go)- This ansatz has a (ferromagnetic) 
solution in the whole T > half-plane, with qo = and zero magnetization for 
H = 0, and qo > and nonzero magnetization for H ^ 0. At zero magnetic 
field and T < 1, this ansatz has another solution with qo ^ that should in 
principle be selected since it has a higher free energy than the first solution. 
However detailed investigation [28j of the eigenvalues of the Hessian shows 
that neither solution is a maximum of the integrand in (|16[) . and accordingly 
both should be rejected, when the absolute value of H is below the so called 
AT line that starts at some nonzero at T = and reaches the H ~ axis 
at Tc = 1 (see Fig[T]). In the rest of the half T > plane, the paramagnetic 
RS solution is stable, which means heuristically that it is the correct solution. 



3 
2 
1 

X 
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-2 
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0.2 0.4 0.5 0.8 1 1.2 

T 

Fig. 1. Phase diagram of the SK model in the T ~ H plane. The paramagnetic RS 
solution is valid outside of the area delimited by the AT lines and the H axis. The 
complementary region is the spin glass phase, where the oo-RSB ansatz gives the 
correct result. 

One is thus led to a more general, non replica symmetric, ansatz. The 
correct (cxd-RSB) ansatz was found by Parisi [11 [7j. It is a hierarchical solution 
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that goes in an infinite number of steps. The first step is the so called "one 
step RSB ansatz". The n replica are partitioned in n/xi blocks of xi replica. 
At this level, Xi is an integer and Xi divides n. One assumes that qa^t = Qi if 
a and b are in the same block, and that qa^t = qo if they do not. With this rule 
we have —§(1 — xi) elements qa.p (for a > b) with value qi, and —^{xi — n) 
elements with value qo- If we plug this information into (j29p . and do formally 
the limit rt ^ 0, this gives 

P{q) = XiS{q - qo) + (1 - Xi)Siq - qi) , (37) 

which means that P{q) has two peaks of weights xi and 1 — Xi respectively. 
Clearly this result only makes sense if xi £ [0, 1]. We started with a, n x n 
matrix subdivided in xi x xi blocks, and then let ri — s- 0. We have to admit 
that in this strange limit, the integer xi became a real variable £ [0, 1]. 

In the next step, the so called "two steps RSB ansatz" , each diagonal block 
of size xi is divided in Xi/x2 sub blocks of size X2 and one assumes that now 
Qa,p = Q2 in the diagonal sub-blocks, and is unchanged elsewhere. 

Performing this procedure k times, we introduce a subdivision n — 
Xo > xi > X2 > ■ ■ ■ > Xk+i = 1, with integer xq, Xi, . . . , Xk, xq/xi, 
Xi/x2, . . . , Xk-i/xk and to a set of values for the overlaps qo, qi, . . . , q^. Then 

k 

P(g)-^(x,+i-x,)%-g.)- (38) 

1=0 

In order for this formula to makes sense, we have to assume that in the 
n — > limit, the Xi variables became = ccg < xi < a;2 < • • • < Xk+i = 1. 
Assuming now that the qiS form an increasing sequence, one has 

/ P{q)dq = ^Y^{xi+i~ Xi) = Xrn+1 m = 0,1, . . . ,k . (39) 

In the k —* oo limit, assuming that the increasing sequence q^'s converges 
to a continuous function q{x), with an inverse x(q), this equation becomes 

Z"'? dx 

/ P{q')dq' = x{q), P{q) ^ ^ ■ (40) 
Jo dq 

We have similarly the equation, for any integer r, 

2 /•! 

y2 9^,6 = yi(-^«+i - X'dli ~ / dxq(x)'' . (41) 

^t^t t^o 

Finally A{q{x)) is expressed as a (complicate) functional of q{x), that 
should be maximized with respect to variations of the function q{x) in order 
to give the free energy /(T): 
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where 0(0, y) is the solution, evaluated at a; = 0, of the equation 



(42) 



a,</>(x,y) = [d^ycl>{x,y)+(3x{dycl,{x,y)f] , (43) 

with the boundary condition 

0(l,y) = /3"Mog(2cosh/3?;) . (44) 

One can show [55] that the (oo-RSB) solution is stable in the whole region 
where the paramagnetic RS solution is unstable . This means (heuristically) 
that we have found the correct solution for all H and T > 0. There is however 
no known close form solution neither for q(a:)r^nor for the corresponding free 
energy f{T), and most articles in the literature use approximate estimates of 
g(x), usually based on the truncated Hamiltonian (see next section), valid for 
T K, Tc- Precise estimates can be obtained however either by an expansion 
of the functional in powers of — T, or by numerical methods (see |30| for 
recent very precise estimates in the H ~ case, and [2J in the H Q case). 

For H ^ 0, the continuous non decreasing function q{x) behaves as follows 
(see Fig[2]): There exist values < Xmin < Xmax < 1 such that q{x) = qmin > 
for X e [0, Xmm], q{x) increases for x e [x,nin, Xmax], and q{x) = qmax for 
X € [xmaxT^- Accordingly the function P{q) has two delta function peaks 
located at q{x) = gmm and q{x) = qmax respectively, and is nonzero between. 
The transition on the AT line is continuous, with qmin ^ Qmax, as the AT 
lined is approached. 

The value Qmax is interpreted as qsA, the overlap between two configura- 
tions in the same pure state 

QEA ^ j;^"^ < '^i >a< >a ■ (45) 

i 

One can show that qsA is independent of a. This definition is in agreement 
with the one given in the introduction ([2]) thanks to the wide separation of 
time scales in the model: A given configuration stays a very long time in the 
same pure state, and accordingly the correlation function C{t) has a plateau 
of height qsA ■ 

As H ^ 0, qmin ^ and the corresponding peak disappears. When H 
is exactly zero, P{q) which is nonzero for < qmin < <? < Qmax for ^ 
becomes symmetric (obviously the two limits N ^ oo and H do not 
commute). At H = 0, P(q) is made of two delta peaks located at ±qmax — 
^Qea with a continuum between. As shown in |30| . and in contradiction with 
many drawings in the literature (and Fig[3]), P{q) has a minimum at a nonzero 



From now on, q{x) denotes the solution of the saddle point equations q'^^{x). 
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value of q (at least as soon as T < 0.96 . . 
with c ^ for small q. When T T^, 
continuous. 



) and behaves like a+bq'^ + c\q\^^ + ■ ■ ■ , 
qEA ^ 0, namely the transition is 



Qmin 



P(q) 



Qmin QEA <i 



Fig. 2. Schematic representation of the oo-RSB solution for H ^ 0: q{x) as a 
function of x and P[q) as a function of q, with qEA ~ qmax- In the actual solution 
the increasing portion of q{x) and the corresponding continuum in P{q) are not 
straight. 



qEA 



P(q) 



qEA 



Fig. 3. Schematic representation of the cxd-RSB solution for H — 0: q{x) as a 
function of x and P{q) as a function of q, with qEA ~ qmax- In the actual solution 
the increasing portion of q{x) and the corresponding continuum in P{q) are not 
straight. In this H — case, P{q) is symmetric in q. Only the q > part of P{q) is 
represented here. 



2.3 Properties of the solution 



Among the fascinating features of Parisi solution, is ultrametricity [32j . Taking 
three clones, one can show that 
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N N 

P(9l,2, 92,3, 93,1) = < %1,2 -E^«'^'^^^'')'^('?2,3 -E'^P^^'") (^G) 

1=1 i=l 

N 

%3,i-E-f^-r')> 

^ ~? TU ^ XI '5((7i,2 - ga.b) 

n(n — 1) n — 2 

^(92,3 - (ib.c)^{q?,,\ - qc,a) ■ 

The 00-RSB solution for the qa,bS has the property that this probabil- 
ity is zero unless the three overlaps 91,27 92,37 93,1 satisfy the ultrametricity 
constraints: at least two overlaps are equal, and the third is smaller or equal 
to their common value. Ultrametricity generalizes to pure states, where one 
defines a distance between pure states a and /3, as 

i 

— Qa.a + QpjJ — 2Qa,^ — 2{qEA — Qa,p) , 

since all self overlaps are equal to qea- 

Another feature [32] of the Parisi solution is that some intensive ob- 

servables are not self averaging El namely some thermodynamic averages 

< A > retain a (disorder) sample dependence even in the large N limit, 
2 

i.e. lim7v-»oo(< A — < A > ) ^ 0. Among self averaging quantities are the 
internal energy, the free energy, the magnetization and the Edwards-Anderson 
order parameter qea- The order parameter distribution function Pj{q) is an 
example of a non self averaging, quantity. The general rule is that observ- 
ables that do not involve correlations between different pure states are self 
averaging, those that involve such correlations are not self averaging. 

In order to show that the order parameter distribution function Pj{q) is 
not self averaging, one considers four clones, and the probabihty ((71,2 7 93,4)- 
Clearly P7((7i.2, (73,4) = Pj[qi.2)Pj{q3A), since clones are non interacting. 
However, we will show that 



^^,7(91,2)^/(93,4) ^ Pj(9i,2) P/(93,4) 7 (48) 

which implies that Pj{q) is not self averaging, since a self averaging Pj{q) 
would not depend on J for a large system. The evaluation of Pj{qi^2, 93,4) is 
done by considering moments of the distribution 



For a discussion of non-self averaging in other random systems, see [33j . 
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J dqi,2 j dqsA 93,4-Pj(9l,2)/'j(g3,4) = < 9l,29|,4 > (49) 



n(n- l)(n-2)(n-3) 



Er s 
qa,bqc,d I 



all different 



with integer r and s. 

We now use the very powerful property of repUca equivalence, namely the 
observation that in Parisi solution one has: 

qa^b independent of a . (50) 

6 

In order to use this property, one notes that . 



ql,bqld = "fa^bqld - 4 Z ia,bql,d (si) 

a,b,c,d a,b c,d a,b,d 
all different 

+ ^Y.1a:bCb ■ 

a.b 

Since the sum b d^a b^a d unrestricted, one can write it as 

E lib E E E <d ■ (52) 



n 

,b d a.b 



In the n limit, one obtains 



1 



< 9^,291.4 > = 3< > + 3< '?l2 > < 9|,4 > : (53) 



Pj{qi,2)Pj{q3A) = ^Pj{qi,2)5{qi,2 - qsA) + 3^^7(91,2) Pj(<Z3,4) , (54) 



in agreement with (|48p . It is interesting that the above relation can be proven 
rigorously |34| (the proof given in this paper is for the r = s = 2 case). 

The method used by Parisi to solve the SK model is far from mathematical 
rigor, to say the less, and it took time to convince the skeptics that it does 
give the correct result. Very recently however beautifully rigorous results have 
been obtained. One is a rigorous proof that the free energy density of the 
SK model has a well defined limit as ^ 00, a long awaited result that is 
far from trivial for a disordered model with infinite connectivity. Then the 
method was extended to show that the free energy of the Parisi solution is a 
lower bound on the free energy of the model |36j , and finally to prove that it 
is equal to the free energy of the model [37] . More recently a proof appeared 
that the AT line is indeed the boundary of the paramagnetic replica symmetric 
region [38] . 
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2.4 Critical exponents 

At zero magnetic filed, the transition at Tc is of second order and one can 
define critical exponents. The order parameter is the mean overlap < q >, 
and one has 



< g > cx (Tc ~Tf T < Tc with /3 = 1 
= T>Tc, 



(55) 



XSG ^N<q^ > , 

cx (Tc - T)-'^ T>Tc 
~ oo T < Tr . 



with 



7=1, 



The values of the critical exponents arc a = — 1 (the specific heat has a 
cusp at the critical temperature), /? = 1, 77 = (the mean field value), 7 = 1, 
and V = 1/2. Hyper scaling holds provided one uses the value d = 6 of the 
space dimension in the formula. This is related to the fact that = 6 is the 
upper critical dimension of the replica field theory. 



2.5 The truncated model 

We have seen that for H = the high T solution has qa,b = Va, b. It is 
accordingly natural to expand A in powers of qa,b in p^ . Only the terms up 
to the order four are usually kept, as one can show that this gives the correct 
critical behavior and the correct qualitative description of the low T phase. 
This truncated model has been very useful historically. It is also quite useful 
in order to study finite size effects (See Subsec. |4T|) . 

The trace over the spins is straightforward. For a given replica index, the 
trace is equal to zero for an odd power of the replicated spins and to two for 
an even power. Neglecting powers of qa.b higher than four, one obtains the 
truncated model (written for H = 0): 



PL 



N 



0^ 1 
In 2 lim In 

4 n^o nN 



n 

.a<b 



I N 
27r/32 



dqab 



exp{-Nf3C[q]) . (56) 



J " ^ XI ^''^^ ^^■'^ 

a.b a,b,c 

~\ ^".^ ^&,c qc,d qd,a + \Y1 ^Ib fa.c " ^ XI '^t'' + ' ' ' ' 

a^b^c.d a.b,c a^b 

where t = (T^ — l)/2. We use the notation q = qP"^ = q/T"^ in order to 
simplify the formula, and write unrestricted sums over the replica indexes, 
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considering {qa,b} as a symmetric matrix witli zero elements on the diagonal. 
For T > Tc, T > and the paramagnetic solution is stable as it should. 

A further simplification, that leads to simpler analytical computations, 
has been proposed by Parisi in one of his seminal articles [5] . It amounts to 
replace ((57)) by 

=^yi ^Iyi ^''■^ ^^-^ ~ 8 ^ ' ^^^'^ 

a.b a,b,c a,b 

keeping the only fourth order term that is responsible for the replica symme- 
try breaking and using the value y = 2/3. This is the Parisi approximation of 
the truncated model, sometimes called the truncated model, sometimes the 
reduced model. In most cases this approximation has only mild effects, chang- 
ing the numerical value of some coefficients. In a few cases however, it gives 
qualitatively wrong results. 

2.6 Some variant of the model 

Two variants of the Sherrington-Kirkpatrick model are worth mentioning. 
The first is the spherical Sherrington-Kirkpatrick model, defined by the same 
Hamiltonian ([1]) as the original model, but with continuous spins ai obeying 
the spherical constraint 

^a.f=iV. (59) 

i 

It can be exactly solved |39j in the thermodynamic limit without the use 
of replica. It can also be solved using the replica trick, and both methods give 
the same results. The physics of the spherical model is however quite different 
from the one of the original model: It is paramagnetic for all T > and H ^ 
but for the region H = and T < 1, where the RS ansatz with nonzero qq 
gives the correct solution. On the other hand, much has been learned of the 
dynamics of spin glass models from analytical studies of the spherical SK 
model (see e.g. [23]). 

The second variant is the p-spin model [40J , where all combinations of p 
spins are coupled together, with the Hamiltonian 

N 

H = - ^ J,i,...^jp cr,i • • -CTip - iJ^cr, , (60) 

l<il<i2<---<4p<JV 1=1 

where the Ji^....^i 's are quenched independent identically distributed random 

p\ 

variables with zero mean and square deviation = 2Np^^ ' "^^^ spins are 
Ising spins cr,; = ±1. The original SK model is recovered forp = 2. Thep > 3 p- 
spin model has a very different physics than the SK model: in the H = case. 
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by decreasing the temperature from T ~ oo one encounters three successive 
phase transitions, a purely dynamical transition, then a transition to a state 
with one-step RSB, and finally a transition to a state with oo-RSB. For details 
seen |40[ I41j . One can define a spherical p-spin model with continuous spin 
variables obeying the spherical constraint [42j [43]. Many analytical results 
have been obtained for the later model. Some subtle qualitative differences 
are found however between the original and spherical p-spin models |41) . 

3 Simulations techniques 

Spin glass simulations are very difficult, since on the one hand the dynamics 
is very slow, and on the other hand one needs to perform the simulation for 
many disorder samples. Both effects are stronger and stronger as T decreases 
and/or N increases. The SK model is no exception. It is furthermore much 
harder to simulate than finite dimension spin glass models since one needs 
0{N) operations to update one spin variable, and not 0(1). The p-spins 
model is even harder with 0{Np~^) operations to update a single spin. For a 
simulation of the p = 3 p-spin model, see [44] . 

The best existing algorithm is called under various names such as "Replica 
Monte Carlo", "Exchange Monte Carlo" and "Parallel Tempering" [15]. This 
algorithm is well known and consists in simulating npx clones of the system 
in parallel at temperatures Ti < T2 < ■ ■ ■ < Tnp^ respectively. Two kind of 
Monte Carlo moves are performed 

• Step 1 (the parallel Metropolis step): Metropolis update of each clone 
independently. 

• Step 2 (the exchange step): Conditionally exchange of the spin configura- 
tions of all pairs of clones, with a suitable acceptance probability. 

In both steps the acceptance probability is the usual Metropolis acceptance 
probability. If tt^. is the Boltzmann weight of state x, with energy E^., namely 
TTj; = exp{—Ex/T)/J2y exp(— ii'y/T), SLadpx,y the probability to go from state 
X at time t to state y at time t + 1, one desires to fulfill the condition 



This condition is called balance (stationarity in the mathematical litera- 
ture). Together with ergodicity (irreducibility in the mathematical literature) 
and the more technical condition of aperiodicity, this ensures that the Markov 
chains converges towards the Boltzmann distribution (351 [JTj |35] . 

Balance is obtained by requiring that the probability to propose the move 
p^Jj and the probability to accept it a^^y satisfy!^ (One has clearly p^.y — 




(61) 



X 




I consider the general case p[ 



7^ Py,x, even if the equality usually holds. 
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^^gp^ ^ min(l,cxp(-l/r(£;. - Ey))) . (62) 

A common misconception is that after the exchange step the system is 
not at equiUbrium, and that this introduces a bias, that one may minimize 
by making many parallel Metropolis steps (step 1) between each exchange 
step. This is not correct, balance is enough to ensure convergence towards the 
Boltzmann weight. 

There is a considerable freedom in implementing this algorithm. The par- 
allel Metropolis step is usually done by updating all spins, either systemat- 
ically (all spins are updated one after the other), or randomly (choose one 
spin at random, then update it, then choose another, . . . ). The exchange step 
is usually restricted to the npT ~ 1 pairs of clones with neighboring tem- 
peratures (the acceptance rate of this exchange is higher), and can be done 
systematically (update the pair at temperatures T\ and then the pair at 
temperatures T2 and Ts, . . . ) or randomly. Step two clearly takes no time 
as compared to step one. One usually alternates the two steps, one parallel 
tempering step (PT step) consists accordingly of one parallel Metropolis step 
followed by one exchange step. 

The choice of the sequence T\ <T2 < ■ ■ ■ < Tnp^ leads also to considerable 
freedom. One constraint is that the acceptance rate of the chain exchange 
(chain of energy E at temperature T with chain of energy E' at temperature 
T') 

=mmil,cxpi-{E/r + E'/T))/exp{-iE/T + E'/r))) (63) 

= min(l - cxp(-(£; - E'){l/T - 1/T'))) , 

should not be too small. The acceptance rate is thus governed by the com- 
bination (T' - TfdE/dT, with dEN/dT = NdeN{T)/dT = A^CAr(T), where 
Cn{T) is the specific heat per spin, which is regular for all temperatures (and 
weakly dependent) in the SK model. One may [JS] adjust for every system 
size N the number npx and the values of the temperature of the clones in 
such a way that (Tj+i — Ti)^7VcAr(T), i = 1, 2, . . . , ?ipT — 1 is roughly indepen- 
dent of i and N . There is however no guaranty that this choice is optimal, for 
example that it minimizes the statistical errors. A more satisfactory prescrip- 
tion is to maximize the number of tunnelings, namely the number of times a 
given clone goes from the minimal temperature to the maximal temperature 
(or the other way), as proposed recently in . 

3.1 An example of Monte Carlo simulation 

In what follows I will use data generated in collaboration with Enzo Mari- 
nari [3] for iJ = (with iV = 64 to 4096, T € [0.4, 1.325]) and Barbara 
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Coluzzi [52] for H ^ (for N up to 3200). For the purpose of these notes, 
I performed additional small simulations in order to measure the number of 
tunnelings as a function of N. I considered N = 64, 256, and 1024. The tem- 
peratures of the clones are equidistant with T = 0.4, 0.425, 0.45, 1.325 
(here npT — 38), with 128 disorder samples (16 for N = 1024), performing 
400000 PT steps, starting from very well equilibrated configurations. 

As shown in Fig. [4l the number of tunnelings is dramatically decreasing 
as N grows. Figure |4] gives another indicator of the algorithm behavior: the 
spread of the distribution of times spent by a given chain at each temperature. 
With infinite statistics, each chain should spend the same amount of time at 
each temperature. I have measured for each disorder sample the maximal (pop- 
max) and minimal (pop-min) time spent (in unit of PT steps) by each chain 
at each temperature. The average time (pop-avg) is equal to 400000/npT = 
10526. These numbers are plotted in Fig. 2] as a function of N. The situation 
is clearly degrading as N grows. Analyzing the N = 4096 data of [Hj (The 
results of a massive numerical effort, with 520000 PT steps and AT = 0.0125) 
one finds pop-min as low as 16.5. 

18000 
16000 
14000 
12000 
10000 

8000 

6000 

4000 

2000 


100 200 300 400 500 600 700 800 900 10001100 

Fig. 4. From top to bottom: maximal (pop-max), average (pop-avg) and minimal 
(pop-min) times spent by a clone at a given temperature as a function of the number 
of sites. Lines have been drawn between the points to guide the eyes. The line in 
the bottom is the average number of tunnelings. 



3.2 Comparison with theoretical expectations 

Figure [S] shows Pj{q) for eight different disorder realizations, using the data 
of [5]. Here TV = 4096 and T = 0.4, a pair of values that is at the borderline 
of what can be done with current algorithms and computer resources. We 
are well inside the spin-glass phase, and the number of spins is equal to the 
number of spins of a 16"^ lattice 0. These eight samples are just the eight first 

This is to say that it is a large system, by current spin glass standards. 
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samples generated by our computer program, they have not been selected 
afterwards. One sees clearly that the shape of Pj{q) is strongly fluctuating 
from sample to sample. In view of (|36p this is very suggestive of the existence 
of several (disorder dependent) pure states. The two extreme peaks correspond 
to the self overlap qea (which is the same for all pure states), the other peaks 
correspond to cross overlaps. If the shape of Pj{q) is strongly fluctuating, the 
position of the outmost peak is not, in agreement with the prediction that qEA 
is self- aver aging. One notes finally that the curves are reasonably symmetric, 
this is a very strong sign that the Monte Carlo statistics is sufficient. 

The non self averaging of Pj{q) makes the measurement of the disorder 
averaged P{q) quite difficult on large lattices. This is exemplified in Fig. |6l 
which presents our estimate of this average with all 256 disorder samples. The 
wiggles are artifacts due to the finite number of samples. 

The Parisi solution has the quite unusual prediction that the low tem- 
perature phase extends to nonzero H (up to the AT line). This prediction is 
unfortunately hard to check numerically, as can be seen in Fig. [7] from [52j . 
Even the very existence of a peak at qmin is not clear from the data. In order 
to see this peak unambiguously, one would need to satisfy two conflicting con- 
strains, H should not be too small, since the weight of the peak goes to zero 
as 7J — *■ 0. It should also not be too close to the AT line since qmin qmax 
in this limit. Clearly much larger systems would bee needed in order to see a 
clear distinct peak at qmin- 

A classical method to locate the critical point from numerical data uses the 
Binder parameter (the Kurtosis of the order parameter distribution function). 
For the SK model this is 



(namely at high enough temperature) and one for a two equal weight delta 
functions distribution. In finite dimension we know from finite size scaling 



that Bn{H,T) is a function of (T - Ta{H))L^/'' = (T - TeXH))N^/^'''^\ and 



accordingly the curves for Bn{H, T) drawn for different values of N should all 
cross at the critical point Tc{H). This provides a very convenient numerical 
method to determine a critical point. In the SK model, finite size scaling holds 
with d = 6 and = 1/2 (see Subsec. \2A\i . The Binder parameter method to 
locate Tc works nicely at zero magnetic field, the curves for different values 
of N decrease monotonously as a function of T, arc well separated away from 
Tj, and cross nicely at Tc- This is not the case at nonzero magnetic field [52j . 
at least with system sizes one can currently simulate. 

It is fair to admit that, from a numerical perspective, the AT line remains 
elusive. This is not a problem in the SK model case since there is no doubt 




When H = < q >= 0, and tlie formula simplifies substantially. 
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5. Pj{q) for eight different disorder realization. Here A'^ = 4096 and T = 
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3.5 



0. 




Fig. 6. The disorder averaged P{q) for q > and zero magnetic field. Here A'^ — 4096 
and T = 0.4. The wiggles are a fluctuation due to the finite (256) number of disorder 
samples. For this value of T, the infinite volume limit of qea is [H] qea = 0.759. 




Fig. 7. The disorder averaged P{q) for q > at nonzero H = 0.3. Here A'^ — 3200 
and T — 0.4. For these values of H and T, the infinite volume limit of qea and q,nin 
are [51] qsA ~ 0.759 and qmin ~ 0.44. The self-overlap qEA is nearly H independent 
at fixed T. This is a consequence of the so called Parisi-Toulouse hypothesis . 

about the existence and exact location of this transition Unc. It become an- 
noying however if one has the EAI case in mind. 

I will close this section by mentioning that there are some numerical evi- 
dences for ultrametricity, in the sense that triplets of typical spin configura- 
tions (with the same disorder configuration) {cr-^^},{crp^},{crp^} do fulfill [53] 
the ultrametric constraint. There is however no numerical evidences |54| for 
the full treelike structure of states, as predicted by Parisi solution. 
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4 Finite size effects for the free energy and the internal 
energy 

Numerical simulations arc obviously limited to finite systems, and simulations 
of spin glasses are indeed limited to very small systems. A detailed understand- 
ing of finite size effects in spin glass models is accordingly highly desirable. In 
what follows, I will only consider the free energy and the internal energy at 
zero magnetic field. I define the exponents fj. and lu according to the equations 



fNiT) = f{T) + B{T)N-^ 
eN(T) ^ e{T) + C{T)N-'^ - 



(65) 



where /(T) and e(T) are the infinite volume free energy and internal energy 
at temperature T. 

4.1 Paramagnetic phase 

In the high temperature phase, the finite size effects for the internal energy 
and free energy can be obtained from the truncated model of (|56l57p . At first 
order, one keeps only the quadratic term in with the result 



/3/^(r) = -ln2-^- lim ^^ln 
4 ri^o nJSI 



0^ . 1 
In 2 lim In 

4 n^o nN 

B"^ 1 
ln2 - ln2T/32 

A AN ^ 



n 

.a<b 



I N 
27r/32 



dqab 



a<b 



n 

.a<b 



2t/32 



exp(-7VT^g2 J 

(66) 
(67) 



The neglected terms in £[q] can be included by expanding the exponent 
as a power series (55l [56l [57] . Each term is represented by a graph without 
external leg. Each line gives a factor 1/{Nt), and each vertex a factor N. 
Introducing the number of lines ^L, of vertices and loops #i? of a given 
graph, one finds that the graph behaves like 7V'^^~^(iVT)~#^ = iV~^^r~#^. 
Organizing the expansion as a loop expansion, one obtains an expansion in 
powers of 1/iV, with the most singular term (as t 0) at each order, given 
by the contribution of the cubic term in C (for which #L = 3(#i? — 1)). The 
development up to the 0{1/N'^) order (four loops) has been obtained in [55] 
(the terms of order are not given in the paper, but have been used in 

the re-summation at Tc). Results up to five loops, for arbitrary n but omitting 
the quartic terms in £ can be found in [57| . 

The internal energy is simply obtained through the equation 
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(68) 



that holds also at finite N. 

As one approaches the critical point both expansions for /jv (T) and e jv (T) 
become singular, and need to be re-summed. Writing symbolically the N 
dependent part of fN{T « Tc) as 



lim In 

n^o nN 



lim In 

n^o nN 



n 

.a<b 



dqab 



exp(-A^(Tg2 + Q3 + Q4)) 



(69) 



.a<b 



2tt 



exp 



(-(Q2+^l/2g3 + ^2/3/^l/3Q4))^ 



where x = 1/(A^t'^), we obtain the N dependent part of fN{Tc) as the x oo 
limit of 



IniV 1 ,. 1 , 
lim — m 

12N N n 



.a<b 



exp(-(g2 + .Ti/2g3 + x'^/^/N^/^Q^)) 



Treating the order four term as a perturbation one obtain finally 



/Ar(Te) = -ln2-l/4 



ew(Tc) = -- 



InTV /(-I) /(-4/3) 
12Ar AT 



-2/3) 



0(iV- 



(70) 



(71) 



2 ■ A2/3 

The behavior of the internal energy is the one expected from scaling. At 
a critical point, the singular part of the internal energy for a L'^ system be- 
havesF^ like L^/'^^'^. Using the values v = \/2 and d = Q (see Subsec. [2^ one 
finds that the singular part of the internal energy does behave like N~'^/^ . In 
order to obtain the values of the coefficients /(-i), and e(_2/3)j one needs to 
handle the theory in the strong coupling a; — > oo regime. This is not easy and 
theses values are poorly determined |55| . 



4.2 Zero temperature 

With the discovery of very efficient algorithms for determining the ground 
state of a spin glass finite system (at fixed J7), the physics of zero temperature 
spin glass has blossomed in the recent years. This includes detailed studies 
of the finite size effects of the internal energy. In the specific case of the 
SK model, good evidences have been obtained for a l/A^^/s behavior for the 
internal energy for both Gaussian and binary disorder distributions |581 1591 
[60l mi [621 [63], as summarized in Tab. [H The value u) = 2/3 is exact for the 
spherical Sherrington-Kirkpatrick model [64j . 

1^ Tliis follows from tfie scaling expression / = 1/L'^F({T - T^)L^'"). 
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P{J) 


N max 




Palassini [59) 


Gaussian 


199 


0.673 ± 0.002 


Bouchaud et al [6U| 


±1 


300 


0.66 ±0.02 


Boettcher [gl] 


±1 


1023 


0.672 ± 0.005 


Katzgraber et al [62] 


Gaussian 


192 


0.64 ±0.01 


Pal [B3] 


both 


2048 


^ 2/3 



Table 1. Numerical estimates for the exponent a; at T = for the Sherrington- 
Kirkpatrick model: reference, disorder distribution, maximum system size and result. 



4.3 Using the Guerra and Toninelli formalism 

As shown in [65j , one can use the so called Guerra and Toninelli interpolation 
formalism, an ingredient in the proof [35] that the free density energy of the 
Sherrington-Kirkpatrick model converges in the infinite volume limit N, as the 
basis of a powerful numerical method to compute the finite size corrections to 
the free energy of the model. Guerra and Toninelli introduced the partition 
function 




(72) 



l<i<j<N 



l<i<j<N/2 



E 

N/2<i<j<N 



exp(/3if^a,) 



that involves a parameter t that interpolates between the SK model with N 
sites {t = 1) and a system of two un-coupled SK models with N/2 sites {t = 0). 
The Jij's, Jj'.j's and J"j's are independent identically distributed Gaussian 
random numbers. It is straightforward to show that 



4 



dt V{t) V{t) < 0, 



(73) 



where 
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1 ^ 2 

i=l 1=1 

(2) ^ 
'^12=]^ 2^ '^'^^ ■ 

i=JV/2+l 

In the above formulas, the di's and r^'s are the spins of two real replica, 
qi2 is the usual overlap, for the N sites system, q[}^ and q^j^2 the overlaps 
restricted to the two subsystems with N/2 sites. The right hand side of ([75]) 
can be evaluated with Monte Carlo simulation. I use the Parallel Tempering 
algorithm, with T e [0.4, 1.3] and uniform AT = 0.025. A total of 2 10^ sweeps 
of the algorithm was used for every disorder sample. The quenched couplings 
have a binary distribution in order to speed up the computer program. Using 
the cumulant expansion ([7]) one can show that in the paramagnetic phase, the 
replacement of Gaussian couplings by binary couplings induces a leading 1 /N 
correction to the free energy density, of the same order as the leading finite size 
correction (see below). One can argue that the effect of the binary couplings 
is also 1/N in the spin glass phase, and is thus negligible with respect to the 
leading l/A^^/s finite size correction: In the spin glass phase the leading finite 
size correction is the same for the binary and Gaussian couplings. Systems 
of sizes N from 128 to 1024 have been simulated with 128 disorder samples 
for each system size (but for N ~ 1024, where I used 196 samples). The 
integration over t was done with the trapezoidal rule, with 39 non uniformly 
spaced points 0. Integrating with only half of the points makes a very small 
effect on the integrand (smaller than the estimated statistical error). The 
data presented at Tc (Fig. [5] and [T^ include the results of an additional 
simulation of a system with N = 2048 sites, limited to the paramagnetic 
phase, with T E [1.0, 1.3], AT = 0.025, with 128 disorder samples, and a 15 
points discretization of t. 

In the low T phase, a remarkable scaling is observed if one plots the ratio 
T>{t)/D{t = 0) as a function of tN'^^^, as shown in Fig. [51 It means that, to a 
good approximation, one has X>(i)/P(0) = F{tN'^^^), with the function F{x) 
decaying faster than 1/x for large x, making the integral in (|73|) converge. 
One has accordingly in the low T phase Jn — foo oc l/N^^^. 

The situation is different at Tc, as shown in Fig.l^l the ratio T){t)/ D{t = 0) 
scales with a different exponent, namely like F{tN^^^), with a large x behavior 
compatible with F{x) oc 1/x, then 

/Ar(r)-/jv/2(r)« / dtV{t)=V{0) f dtF{tN^/'') (75) 
Jo JO 

cx 1/iV IniV/iVo , 

^'^ With a distribution adapted to the shape of 'D{t), that is peaked at low t. 

Admittedly much larger system sizes would be needed in order to be sure that 
the system really obeys this asymptotic behavior. 
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Fig. 8. T>{t) /D{t = 0) as a function of tN^^^ (both in logarithmic scale) , for T = 0.6. 
The full line shows the 1/x behavior, the dotted line shows the behavior. 
Clearly D{t) decreases faster than 1/x for large x. The precise behavior of V{t) is 
not essential for my argument, as soon as T>{t) decays faster than 1/x. 



for some undetermined Nq. Use has been made of the fact that, according 
to finite size scaling, V{Q) ^ -l/2< {q[]l)^ > = -l/2< (qj^^)2 > scales like 
N^0/(d>^)G{{T - Tc)Afi/('*'')), with a finite non zero hm^^oG(a;) and that 
P/{dv) = 2/3 (see SubsecE!]). This behavior of /jv(T) - fN/2{T) is in agree- 
menlEl with (ffO]) . 



128 


X 




256 


□ 




512 







1024 






2048 


V 




2./X - 




V 



0.001 0.01 0.1 1 10 100 

tN'"- 



Fig. 9. V{t)/D{t = 0) as a function of tTV^/^ (both in logarithmic scale), for T = Tc. 
The straight line shows the expected 1/x behavior, in order to guide the eyes. 



Figure [TUl shows, as a function of l/N'^^^, my estimate at T = 0.4 of the 
difference (/at — fN/2)/T, obtained by integrating numerically ((73)) . compared 
to the result of a hnear fit [Jn - fN/2) /T = -A/N-'^/^, with A = 0.82 ± 0.02 



But the value of the prefactor is not predicted by my method. 
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Fig. 10. Numerical data for (/jv — /jv/2)/r as a function of together with a 

numerical fit to the data of tlie form {fN-.fN/2)/T = -A/N^^'-', with A = 0.82±0.02 
dotted line. Here T = 0.4, N = 128, 256, 512 and 1024. 



and — 4.9. The agreement is good within estimated statistical errors. A 
similar agreement is obtained for other values of T in the spin glass phase, 
e.g. A = 0.39 ±0.01 with = 3.6 for T = 0.6, and A = 0.18 ±0.01 with = 
33 (a large value presumably related to the proximity of the critical point) 
for T = 0.8. The value = 2/3 is in contradiction with the old analytical 
prediction of [66| . In this paper, arguments arc given that relate lo to the 
exponent p of the first correction term in the expansion of the replicated free 
energy — 1/ (3lmiN^oc> 1 /N In Zj in powers of n, by the equation uj ~ (j)^l)/p. 
In the paramagnetic phase it is known that there is no term in this expansion 
beyond the linear term, i.e. p is infinite, and thus a; = 1 in agreement with 
the previous discussion. In the spin glass phase Kondor [67j has found, using 
the tnmcated model of ([5^ . that p = 6, and thus a; = 5/6. The resolution of 
this contradiction lays presumably in the use of the Parisi approximation (j58[) 
in [67], and we can conjecture from our data that indeed p = 3. 

Since the energy at zero temperature (and thus the free energy) also be- 
haves like l/A^^/"^, the most natural conclusion is that the leading finite size 
corrections for both fN{T) and eiq{T) behave like in the whole low 

temperature phase, for both binary and Gaussian distributions. This seems 
however to contradict the numerical results of [68| . In this paper, the internal 
energy of the SK model with Gaussian P{J) has been measured with Monte 
Carlo simulations for values of T between T = 0.1 (a very low value made 
possible by the small sizes simulated) and T = 1.22, with N = 36, 64, . . . , 196. 
The data for eAr(T) are fitted as eAr(T) = e{T) + a{T)N-'^, with a result 
for UJ that is compatible with 2/3 for both T = and T = Tc, but with a pro- 
nounced deep between, in contradiction with my conjecture. Using the data 
for eAr(T) obtained from the simulation of [9], it is quite simple to repeat the 



Using the precise estimates of e{T) from |30) . 
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analysis of [55] for systems up to iV = 4096 (but with binary P{J)). Figure [TT] 
shows my results compared to the one of [55] ■ The use of larger system sizes 
clearly reduce the effect observed in this paper. The most likely conclusion is 
that the analysis of [SQ is affected by systematic errors due to the small sizes 
used. The shape of the sub-leading corrections to the energy is not known 
below Tc, it is known at the critical point, however, and there the sub-leading 
corrections are decaying very slowly, the expansion is indeed an expansion in 
powers of and it is may be not so surprising to have difficulties finding 

the correct leading exponent from data with 36 < iV < 196. 

1 



0.9 
0.8 

3 

0.7 



0.6 
0.5 

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 
T 

Fig. 11. Behavior of the finite size exponent uj as a. function of temperature for the 
Sherrington-Kirkpatrick model. From top to bottom: results of a (gnuplot) fit of 
the data of [5] with 4096 > N > 256 (with estimated errors); fit of the same data 
with 128 > TV > 64 (without estimated errors, in order not to clutter the figure) 
and data from |68] . 

Figure [T2l shows my estimates for {/n — fN/2)/T at Tc as a function of 
1/A^, together with the prediction of f7D|) . A good agreement (with = 4.3 
if one excludes the N ~ 128 data from the fit) is obtained using the value 
= 7.8 ± 0.2, namely = ln(7.8)/12 = 0.17. . ., within the range of 

results presented in [55] , 

The method has been extended [69] to the computation of the fiuctuations 
of the free energy, with the result that, in the spin glass phase, A"^ f^iT) = 

f^{T) - fN{Tf oc N-"^", with a w 3/5. This is definitively smaller than the 
value found atr = 0[59l|60l|6ll[6l[S3], that is a ^ 3/4. 

4.4 Conclusions 

The picture that emerges is the following: above Tc the finite size corrections 
of both fN{T) and e]\f{T) are given by series in powers of l/N, whose terms 
can be evaluated by perturbation theory. This expansion becomes singular 
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Fig. 12. Numerical data for (/jv — /jv/2)/r as a function of 1/N, togetiier witii 
the beiiavior implied by tlie equation: /iv/T = foo/T + l/(12iV) XuN/Nq. The full 
line is drawn with the value A'o = 1. The dotted line is drawn with the value 1/7.8, 
from a fit to the data. Here T = 1, N = 128, 256, . . . , 2048. 

as the critical point is approached. At the critical point the size dependent 
part of the free energy behaves hke (l/12iV) InA^/A^o and the size dependent 
part of the internal energy as e(_2/3)/A^^^'^, with coefficients that are poorly 
determined analytically, but can be evaluated with Monte Carlo simulation. 
The sub-leading corrections are decaying very slowly with N . In the spin glass 
phase, the data are compatible with a l/A^^/"^ behavior for both the leading 
finite size corrections to e^iT) and /jv(r). 
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